General heatbath algorithm for pure lattice gauge theory 
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A heatbath algorithm is proposed for pure SU(JV) lattice gauge theory based on the Manton 
action of the plaquette element for general gauge group N. Comparison is made to the Metropolis 
thermalization algorithm using both the Wilson and Manton actions. The heatbath algorithm 
is found to outperform the Metropolis algorithm in both execution speed and decorrelation rate. 
Results, mostly in D = 3, for N — 2 through 5 at several values for the inverse coupling are 
presented. 
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I. INTRODUCTION 

A heatbath algorithm is proposed for pure SU(iV) lat- 
tice gauge theory based on the Manton action of the pla- 
quette element [l| . The relation of the Wilson and Man- 
ton actions is equivalent to the relation of the Frobenius 
and Riemann metrics for the distance between elements 
of the gauge group @. While heatbath efficiency has 
been achieved for the cases of SU(2) and U(l) using a 
biased- Metropolis algorithm 3|, 5| > an d a heatbath algo- 
rithm for SU(2) is known [H,[6| which may be extended to 
N > 2 using covering subgroups Q , the composition of a 
direct heatbath updating scheme for general gauge group 
N has been an outstanding problem in lattice gauge the- 
ory for quite some time. As the direct approach seems 
intractable, we consider an indirect means to accomplish 
the thermal updating. By relating the Riemann norm of 
a plaquette in the matrix representation Q to the Frobe- 
nius norm of the plaquette in the vector representation q, 
one may generate random plaquette elements of known 
Manton action, which reduces to the Yang-Mills action 
in the continuum limit. The algorithm here generates 
a proposal for the updated mean plaquette by averag- 
ing some number of random elements with action drawn 
from a gamma distribution. Dividing out the mean en- 
vironment contribution then provides the updated link 
value. 

After calibration against known results [8| , we compare 
the performance of the proposed heatbath algorithm to 
the Metropolis thermalization algorithm Q using both 
the Wilson and Manton actions [lfj. We show that one 
may generate identical lattice action distributions using 
either thermal updating scheme, with demonstration for 
a D = 3 and L = 4 isotropic lattice having periodic 
boundary conditions. The running time of the heatbath 
algorithm compares favorably with that of the Metropolis 
algorithm using either action and displays the expected 
linear dependence on lattice volume L D . Both the link 
and plaquette decorrelation rates are evaluated for a par- 
ticular case of N and f3, where the heatbath algorithm 



is found to outperform the Metropolis algorithm in both 
execution speed and decorrelation rate. 

Summary of notation: 
Isotropic Euclidean lattice regularization i/a6 {L D } 
Special unitary group SU(iV) with identity Tr I = N 
Link variables U £ {U}, plaquette variables Q £ {Q} 
Exponential parametrization Q — exp(i ^ a q a H a ) 
Generator normalization Tr H a Hb = 25 a b 
Conjugate transposition = U~ x , = H 



If 



Frobenius norm ||Q|| 

Frobenius metric dpiQiiQz) — WQi — Q2I 
Riemann norm = ||logQ|||,/2 

Riemann metric d R (Qi,Q 2 ) = ||log(QJ Q2)\\f/V% 
Boltzmann factor exp(— S) = exp(— PY^q Sq) 
Wilson action S W {Q) = d|(7, Q)/2A^ = 1 - Re Tr Q/N 
Manton action S M (Q) = d 2 R (I,Q)/N = ||logQ||f./2iV 
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II. COMPARISON OF ACTION DEFINITIONS 

The Wilson action is recognized as a measure of the 
gauge invariant plaquette element's distance from the 
identity normalized by the gauge group TV. Working 
backwards from the usual expression, one writes 

Sw(Q) = 1 - Re Tr Q/N , (1) 
= Tr [27- (Q + Q+)]/2Ar, (2) 

= Tr [(/-Q) t (/-Q)] /2N , (3) 

= d%(I,Q)/2N , (4) 

where the Frobenius metric measures the distance from 
the identity to the plaquette in the space of the gen- 
eral linear group GL(N, C). Similarly, the Manton ac- 
tion is related to the Riemann metric of the plaquette 
matrix Q, which measures the distance along the one 
parameter subgroup of SU(A r ) connecting I to Q, given 
by Q p = exp(plogQ) for < p < 1. Using the expo- 
nential parametrization Q = cxp(i^2 a q a H a ) for trace- 
less Hermitian generators H normalized to Tr H a Hb = 
28 a b, the Riemann norm in the matrix representation is 
equal to the Frobenius norm in the vector representation, 
lllogQUl = Tr E a <l 2 a H a = 2||?|||, so that the Man- 
ton action is proportional to squared radius of the pa- 
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FIG. 1. Displayed as + is the ratio of the matrix and vector 
norms ||Q||ji/||?||f averaged over 100 random unit vectors q 
for various N. Also shown as is the ratio of the Wilson 
and Manton actions Sw(Q)/Sm(Q) which approaches unity 
as \\q\\F —> 0. The injective limit ||g||p = n corresponds to an 
element antipodal to the identity. 



rameter vector q = rq, giving Sm(Q) °c ||<5||jj = r 2 . 
The maximum for either action is given by its value 
for an element antipodal to the identity Sq(I'), where 
I' = diag[— 1, . . . , — 1, ±1] and the final sign ensures 
detl' = 1, yielding S W (I') = 4[(N/2)/N and S M (I') = 
[(N/2)tt 2 /N for [(N/2) the greatest integer < N/2. 

The exponential parametrization of SU(iV) elements 
is unique out to the injective limit r < 7r — beyond that, 
the matrix iY] a q a H a is not the principle logarithm of the 
corresponding plaquette Q, indicating that a shorter path 
through the group connects the element to the identity, as 
seen in Fig.[TJ For reasonable values of /3 close to the con- 
tinuum limit, the chance of encountering elements near 
the injective limit is exponentially suppressed. Writing 
the logarithm as logQ = — X)^Lr(^ — Q) k l^i we see that 
the Wilson action corresponds to taking only the first 
term in the expansion. Consequently, it may be viewed 
as the leading order approximation to the Manton ac- 
tion, similar to the relationship between the chordal and 
arc distance between two points on a circle. In the limit 
of vanishing action Sq — > 0, the two expressions become 
equal, so that the ratio Sw(Q)/Sm(Q) —> 1, also shown 
in Fig. [T] — at finite j3, one expects Sw{Q) < Sm(Q)- 

The equivalence of the action definitions in the contin- 
uum limit may be shown analytically. The relation of the 



Wilson action to the Yang-Mills action is well known [lfj 
and will not be repeated here. Picking up the deriva- 
tion for the Manton action in D = 4 at the plaquette 
discretization Q^ v (x) — exp[— a 2 G^ v (x)} for G fiv (x) — 



F liv (x) + 0(a), one finds logQ 
one may write 



tfF^+Oia?). Then 



Sm(Q) = ^dl(I,Q) = — \\\ogQ\\ 2 F = — Tr G\ V G^ , 

(5) 

and doubling the plaquette sum 2J2q = Sa; ^ ^ one 
tains 



Q 



S = PY, S m{Q) = -^£a 4 iV(z)i^>) + 0(a 5 ) 



(6) 

where the minus sign appears because A^(x) — —igA b ^Hb 
for coupling constant g. Thus, pure lattice gauge field 
theory may be described equally in the limit j3 — > oo by 
either the Wilson or the Manton action. The coupling 
may become dimensionful (3 = 2N/a 4 ~ D g 2 when D ^ 
4. The numerator for f3 is conventional and serves to 
cancel the normalization factor in the denominator of 
the definition of the action; it might be more physical 
to normalize by the number of degrees of freedom in the 
gauge field d = N 2 — 1 so that the normalized Manton 
action becomes the mean squared value of the parameters 
in the vector representation. 

Later we will make use of the average of some number 
of SU(JV) elements, so let us now look at how the pro- 
jection of the sum relates to the choice of action. The 
link to be updated U is surrounded by uq = 2 D ~ 1 sta- 
ples Gk € SU(A r ) composed of the remainder of the pla- 
quettes involving U whose sum belongs to the general 
linear group ^2 k Gk — G £ GL(N, C). Generalizing the 
discussion by Moakher [2j, the Frobenius mean is the 
group element Gf £ SU(iV) which minimizes the met- 
ric of GL(AT,C), such that G F = argminy^. d 2 F (U, G k ), 
which one may show is equal to the projection of G onto 

SU(iV) that maximizes Re Tr G F G. The most expedient 
evaluation [lj is first to project from GL(N, C) to 
SL(N, C) by dividing out the A^th root of the determi- 
nant G = G/IGf/^and then from Sh{N,C) to SU(A^) 
using Gf = G(G^G) -1 / 2 . Alternately, one may write 
G F = G(G^G)- l / 2 \G- 1 G^\ 1 / 2N , requiring care that the 
determinant comes out 1 when G is far from SU(AT). The 
Riemann mean is the element Gfi which minimizes the 
metric of SU(A r ), such that Gr = argmin^ fc d 2 R (U, Gk), 
which coincides with Gf for the case uq — 2; otherwise, 
the means may differ but only by some small amount 
for elements sufficiently close together. No simple for- 
mula exists for Gr in the general case (to our knowledge), 
though it may be found from a numerical minimization 
of the metric, and so for the following we will write G 
for the projection of the staple sum, using Gf for its 
evaluation when we really ought to be using Gr. 

The bi-invariance of either metric ensures that dis- 
tances are preserved under a unitary transformation 
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FIG. 2. Distances are preserved under a unitary transformation of the group elements. The environment projection G minimizes 
the sum of the distances to the elements Gk, as shown by the dashed lines. The distances appearing in the action are shown 
by the solid lines. 



d(UAV,UBV) = d(A,B), which for the Wilson action 
accounts for the permutation of the order of the elements 
under the trace. Thus, one can relate graphically the ac- 
tion of the plaquettes Q k = UGk to the distances from W 
to the environment elements Gfc, as shown in Fig. [2] The 
projection of the plaquettes is given simply by Q = UG. 
One can see that the action of the projected plaquette 
sum need not be simply related to the sum of the pla- 
quette actions. 



III. DESCRIPTION OF THERM ALIZATION 
ALGORITHMS 

We will start by discussing the commonly used al- 
gorithm by Metropolis as well as the Creutz-Kenncdy- 
Pcndlcton heatbath algorithm for SU(2) and why its di- 
rect extension appears intractable for N > 2. We will 
then describe our prescription for a general heatbath al- 
gorithm based on the Manton action. We will show in the 
following section that identical action distributions are 
generated by the Metropolis and heatbath algorithms, 
with some allowance for the case of N = 2. On a techni- 
cal level, we have found that, even for the case of SU(2), 
while the overrelaxation algorithm [l3T-[l5l] leaves the Wil- 
son action invariant for U — > G^pU^G^p, the Manton ac- 
tion is not invariant (though is nearly so for G R ), and so 
we will not be including microcanonical overrelaxation in 
this study. 



A. Metropolis algorithm 
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The usual implementation of the multi-hit Metropolis 
(MP) algorithm updates a link U ncw «— U a id by accu- 
mulating proposals P <G SU(iV) on the left of the original 
link U DCW — PU id- If the proposals P are chosen such 
that inverses are equally likely pt{P) — Pt(P ), then 
the trial probabilities cancel in Eq. (O, leaving an ac- 
ceptance probability proportional to the Boltzmann ra- 
tio pa oc exp[-/3(S' ncw - S id)] = cxp(-/3As), which 
is compared to a uniform deviate p G [0,1]. Here, the 
hp = 10 proposals P = exp(z ^2 a p a Ha) are generated by 
a Gaussian deviate on the parameters p in d dimensions 
with zero mean and variance controlled by a parameter 
e p <C 1, which adjusts during thermalization to tune the 
hit ratio Ha/ hp ~ 1/2 then is held fixed during measure- 
ments, so that P and P' have equal likelihood peaked at 
the identity. 

Using the Wilson action, one can write Acf = 
Re Tr (I - P)U oid G/N, where G is the sum of the sta- 
ple elements Gk surrounding the link to be updated. 
For the Manton action, each term in X)fc'^Af(Qfe) must 
be evaluated; the most efficient means is to notice that 
Sm(Q) = — X)j(l°g^i) 2 /2^ where Xj are the eigenval- 
ues of Q = VqAqVq. The additional calculational load 
approximately doubles the running time of the MP al- 
gorithm using the Manton action in D = 3. We verify 
our construction of the MP algorithm by comparing its 
results to those presented by Teper [al making use of the 
algorithm by Cabibbo and Marinari [7|. 



The benchmark algorithm for the thermal updating of 
gauge field configurations according to the Boltzmann 
factor exp(— S) = exp(— PJ2q &q) = II q Wq is that 
of Metropolis @, where the local transition probabil- 
ity for non-uniform sampling is composed of a trial and 
acceptance probability p(U ncw f7 ld) = PA(U ncw 



B. Creutz-Kennedy-Pendleton algorithm 

For the case of N = 2 a direct heatbath algorithm 
is known [5, 6]. The CKP algorithm makes use of two 
features peculiar to SU(2): the simplicity of the spectral 
parametrization [161 ] of group elements, U = uqI + iu- H 
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FIG. 3. Generalization of the CKP algorithm in terms of the 
group metric. The proposal P r is parametrized by its distance 
from the element which minimizes the action. These elements 
are related to those of Fig. [2] by Gk = G (G^Gk)- 




FIG. 4. Indirect heatbath updating generates ug proposals 
for plaquettes according to the canonical distribution whose 
projection is P. Supposing Q — > P given fixed G requires the 
updated link to take the value (/ nC w = PG . 



for H = 3 the Pauli matrices and uq — ±(1 — ||u|||') 1 ^ 2 , 
and the proportionality of a sum of elements to another 
group element, G = G^G] 1 ' 2 . The Haar measure is re- 
duced to dll = duo(l — u 2 ] ) 1 / 2 d 2 Q u , where Sl u is the solid 
angle along u, and the Boltzmann weight according to 
the Wilson action is Uq w Q « exp(/3Tr UG/2). The 
action is parametrized as UG = {G^UG = \G\ 1/2 Q 
for Q = qol + iq ■ a, yielding a final distribution 
JdQU Q W Q cx /^^(l-^^exp^lGl^go)/^. 
The angular distribution may be generated from normal- 
ized Gaussian deviates as above, leaving the single pa- 
rameter q to be determined. The Creutz algorithm gen- 
erates the exponential factor from the cumulative distri- 
bution and corrects for the Haar measure using a rejec- 
tion step, while the Kennedy-Pendleton variant includes 
a piece of the Haar measure with the exponential distri- 
bution under a suitable change of variable and corrects 
for the remainder with a rejection step. 

The generalization of the CKP algorithm is depicted 
in Fig. |3J related to the right hand side of Fig. [5] by 
Gfc = G(G^Gfc). The proposal is parametrized as P r = 
exp(zrp ■ H) such that P Q = I and P~ r = P r \ where 
p is a random unit vector and r is the Riemann met- 
ric between that element which minimizes the action Iq 
and the proposal P r . Decomposing the angular depen- 
dence P 1 = V\ Ai Vi lets one write the sum of the actions 

as S sum — J2k Sq(ViA^ t Vi G^ Gk), which one relates to 
an action value drawn from the desired distribution by 
Sdraw — Ssum /no- Using the Wilson action, one can 
write =flG -Re Tr P-^G^G) 1 ' 2 /N\G~ l G^\ l l 2N 

where (G^G) 1 / 2 is not an element of SU(AT). For the 
Manton action, the evaluation of the sum as a function of 
r is even less tractable, so that no clear means of generat- 
ing its distribution presents itself. We remark, however, 
that a biased-Metropolis algorithm 0, 0] may be built 
along these lines, where the distance r would appear in 
the trial probability pt{P t ) and the distances for S sum 
appear in the Boltzmann ratio. 



C. Indirect heatbath algorithm for general N 

The first step in developing a local heatbath (HB) algo- 
rithm is to break the relationship between the links U new 
and [/ id so that updates are generated with the canon- 
ical distribution dU aew exp(— /3S new ) without regard to 
the current link value As the direct generalization 
of the CKP heatbath algorithm is intractable, we pro- 
pose an indirect scheme for generating the updates U ncw . 
We start by generating uq proposals for plaquette ele- 
ments Pk with action drawn from a gamma distribution 
as detailed below — one may think of these as selecting 
new environment elements for a fixed link value of J7 id- 
After taking their projection P, we reverse the sense of 
what has changed, supposing the environment projection 
remains G while the plaquette projection Q — > P, so that 

the updated link takes the value U ncw = PG , as shown 
in Fig. HI While it would be nice to get P with the 
generation of a single random element, the action of the 
projected sum does not appear to be simply related to 
the sum of the actions for TV > 2. 

The generation of the gamma distribution for the 
action Sm(Pic) is accomplished using the rejection 
method [TtJ • Retaining only the radial coordinate, the 
Haar measure fioj ] (accounting for our generator normal- 
ization) is dP oc dr r d ^ 1 cxp[— r 2 N/6 + C(r 4 )], and the 
local partition function Zp is written as 



,„>, ■■■>„■/'. x drr d-l e -(fi/N+N/^ j (g) 



ll-r'lh 

drr a e- br , 

o 



- (^,b\\ni 



(9) 
(10) 



where j(a,y) = f~ dr r a ~ 1 e~ r is the incomplete gamma 
function. The cumulative distribution is F(r) — "/[(a + 
l)/2, br 2 ]/Zp, and the unnormalized frequency distribu- 
tion is /(?') = r a exp(— br 2 ) oc dF/dr. The frequency 
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FIG. 5. Normalized frequency and cumulative functions for 
the gamma (solid) and Gaussian (dashed) distributions as 
described in the text. Also shown as x is the frequency dis- 
tribution histogram of 10 4 generated values. 



distribution is peaked at r p = (a /2b) 1 / 2 with a value 
f p = (a/26) a / 2 exp(-a/2). The distribution is approxi- 
mated by first fitting a Gaussian to the range r > r p by 
minimizing J dr {exp[— (r — r p ) 2 /2a 2 } — f / f p } 2 and then 
writing g(r) = / p exp[— (r — r p ) 2 / s 2 ], where s = (3d 2 ) 1 / 2 
expands the Gaussian so that g > f. Its integral yields 
the cumulative distribution 



G(r) 



dr'g{r') 



crf(^) -erf Tp 



(11) 
(12) 



which one may invert analytically. For a value Gdraw 
drawn uniformly from the range [0, G(|| /'Hn)], one gets 



rdr 



r p — s erf 1 



erf(^ 

s 



Gdr 



(13) 



which is accepted with a probability of f (rd ieLW ) / g(r^ raw ) . 
The normalized values of these distributions are shown 
in Fig. [S] for N — 3 and f3 = 24, along with the fre- 
quency distribution of 10 4 draws. Due to the close ap- 
proximation g ~ f, there is a fairly high acceptance rate 
of F max / G 

max > 75%. 



IV. LATTICE ACTION DISTRIBUTION 

To evaluate the performance of the indirect heat- 
bath algorithm, we compare the expectation value of 
the plaquette action using the Manton definition Sm = 
(Sm{Q)){q} on a D = 3, L = 4 isotropic lattice measured 
over 100 configurations separated by 10 thermal sweeps 
following 100 thcrmalization sweeps from a warm start. 
We consider the MP algorithm using either the Wilson 
or Manton action, and for the HB algorithm we compare 
the effect of neglecting the exponential factor of the Haar 
measure bo = j3/N with its inclusion bn — bo + N/6, 
as shown in Fig. [HI Normalizing the coupling by the 
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FIG. 6. Comparison of the mean Manton action value Sm for 
the Metropolis and indirect heatbath algorithms on a L D = 
4 3 lattice as a function of the normalized coupling for iV 6 
{2, 3, 4, 5} plotted as {□, 0, O, +}■ The mean action for MP 
using Sw is higher than for MP using Sm, as is the value for 
HB using bo compared to using bn- The SU(2) values for HB 
are consistently higher than expected. 



number of degrees of freedom, the N dependence of the 
mean action disappears so that the values of S , m(/3/^) 
all follow the same linear relation on logarithmic axes — 
a similar plot obtains from the values presented by Te- 
per @ . The effect of the choice of action definition in the 
MP algorithm is most noticeable at low /? values, where 
SM(MPff) > S M (MP M )- For the HB algorithm, we find 
that the mean action for SU(2) is consistently slightly 
higher than the corresponding values for N > 2; we have 
not yet been able to root out the source of the discrep- 
ancy, but as the direct CKP algorithm applies here this 
issue is of no great concern. Interestingly, we find agree- 
ment between using HB with bo and MP with Sw an d 
between using HB with bn and MP with Sm > which is un- 
derstandable considering that both bn and Sm include 
higher order corrections to bo and Sw a s elements get 
further from the identity. 

We next compare the plaquette action distribution for 
SU(3) at P = 24 thermalized by MP with S M to that 
for HB with displayed in Fig. [7] as histograms using 
either linear or logarithmically spaced bins. The relation 
between the linear and logarithmic histograms may be 
expressed [H, E^l as JqASm = SMfQd(\ogSM), and a 
least squares regression over the logarithmic bins of a 
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FIG. 7. Frequency distribution of the Manton plaquette ac- 
tion for SU(3) at /3 = 24 plotted as histograms x with either 
linear or logarithmically spaced bins. The regression func- 
tion Eq. (TT4)) is shown as the solid line. The MP and HB 
algorithms return virtually identical action distributions. 



TABLE I. Coupling values extracted using Eq. (|14|l from pla- 
quette action distributions generated on a L D = 4 3 lattice. 



N 


P 


2 


6.00 


9.00 


12.00 


15.00 


6 


16.00 


24.00 


32.00 


40.00 


4 


30.00 


45.00 


60.00 


75.00 


r. 
O 


48.00 


72.00 


96.00 


120.00 


N 


/Sat for MP with S w 


2 


5.55 


8.50 


11.62 


14.72 


6 


14.48 


22.52 


30.61 


39.12 


4 


26.87 


42.02 


57.47 


72.30 


r, 
o 


42.96 


67.30 


91.56 


116.10 


N 


/3fit for HB with bo 


2 


5.11 


7.99 


11.07 


13.85 


Q 

o 


14.65 


22.78 


31.17 


39.23 


4 


27.38 


42.60 


58.24 


73.19 


r. 
u 


43.44 


67.80 


92.29 


117.10 


N 


/3at for MP with Sm 


2 


6.13 


9.30 


12.14 


15.27 


3 


16.08 


24.42 


32.26 


40.77 


4 


30.28 


45.81 


60.93 


75.75 


5 


48.62 


72.78 


97.22 


121.67 


N 


(Sat for HB with 6^ 


2 


5.70 


8.65 


11.70 


14.52 


3 


16.28 


24.24 


32.57 


40.86 


4 


30.13 


45.36 


60.66 


76.79 


5 


47.90 


72.07 


96.53 


120.42 



RUNNING TIMES AND DECORRELATION 
RATES 



function 

Sq « S M +d ' 2 ex P (-fotS M D/2) (14) 

yields an estimate of the coupling from the generated ac- 
tion distribution of ^(MP^) = 24.42 and /3 fit (HB ff ) = 
24.24. (Technically one should perform a Bayesian re- 
gression to extract the parameters, but for our purposes 
here a simple fit is adequate.) The fitting function also 
is shown in Fig. where the agreement of the action 
distributions is apparent. 

We use the fitting function Eq. (|14|l to extract the 
coupling from the generated distributions used in Fig. [5] 
for N G {2,3,4,5} and (3/d G {2,3,4,5}, displayed in 
Table Q] For consistency of comparison we use Sm in 
Eq. ([T3]) even when using Sw or 6 in the thermaliza- 
tion algorithm — these /3at are lower as Sw < Sm- The 
/3fit values for MP with Sm and HB with bjj are in ex- 
cellent agreement with each other and with the coupling 
/3 used in the updating algorithm, except for the case 
N = 2 which seems to show a slight systematic discrep- 
ancy as mentioned above. Certainly for the larger N, 
the Metropolis and indirect heatbath algorithms gener- 
ate mutually compatible lattice action distributions when 
using the Manton action and accounting for the exponen- 
tial factor of the Haar measure. 



To compare the efficiency of the algorithms, we mea- 
sure the running time t±o of 10 thermalization sweeps for 
lattices with L G {4,6,8} and D G {2,3,4}, as shown 
in Fig. [51 No optimizations have been made beyond the 
obvious, and the evaluation is on a serial processor. The 
expected linear dependence on lattice volume tio cx L D 
is observed. The ratio of the running time for the MP 
algorithm using the Manton action to that for the Wil- 
son action increases along with the lattice dimension D. 
The running time for HB (using bu) is always less than 
that for MP using Sm and is usually less than that for 
MP using Sw- The indirect heatbath algorithm allows 
for the use of the Manton action with a computational 
efficiency exceeding that of the corresponding Metropolis 
algorithm. 

To measure the correlation between two configurations 
{U(x)} and {U'(x)} separated by a given number of 
update sweeps, we consider both the link correlation 
C(U,U') = 1 - \\WU'\\ R /\\r\\ B . and the plaquette cor- 
relation C{Q,Q') = 1 - WQ^Q'Wr/WI'Wr evaluated for a 
particular site and averaged over the lattice directions 
and volume. This measure of correlation equals 1 when 
U = U' and descends to some residual value determined 
by the inverse statistical temperature, as for high /3 none 
of the elements may stray far from the identity — early on 
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FIG. 8. Running times tio for 10 thermal updating sweeps as a function of lattice volume L D for MP with Sm shown as □, for 
MP with Sw as (), and for HB with bu as Q- 



we reproduced values similar to those of Ref. [13J using 
that definition of correlation. 



We evaluate the link and plaquette correlations on a 
L D = 4 3 lattice with f3/d — 3, comparing MP with Sm to 
HB with b H averaged over 100 measurements. Only for 
N = 5 was a slight volume dependence seen for MP com- 
pared to a short run with L = 6. The link correlation 
C(U, U') decays much faster for the HB algorithm, as 
shown in Fig. [9j where the rate of decorrelation displays 
a dependence on the gauge group N. As links are not 
gauge invariant objects, we also measure the plaquette 
correlation C(Q, Q'), shown in Fig. [TUJ Here we see that 
the HB algorithm has achieved nearly the minimal corre- 
lation after only a single sweep, in sharp contrast to the 
MP algorithm which requires a greater number of sweeps 
to reach minimal correlation as N increases. Having bro- 
ken the relationship between the links U ncw and U a \d in 
the thermal updating, the indirect heatbath algorithm 
outperforms the Metropolis algorithm at decorrelating 
the lattice gauge field configurations. 



VI. LINK VS PLAQUETTE SIMULATION 

As an interesting conjoinder to our investigation of 
link updating by the heatbath algorithm, we would like 
to present some thoughts on an alternative approach to 
quantum gauge field simulations on a Euclidean lattice. 
An efficient means of generating random elements Q with 
known action Sm(Q) suggests the possibility of running 
a lattice simulation with plaquettes rather than links as 
the dynamic variables. As there is a 1 : 1 geometri- 
cal correspondence between the number of links and the 
number of plaquettes, a given lattice configuration {U} 
may equally well be represented as a configuration {Q}. 
Furthermore, the representation in {Q} is gauge invari- 
ant, corresponding to an entire class of {U} related by 
gauge transformations. Essentially, the plaquette repre- 
sentation "integrates out" the gauge degrees of freedom 
to leave one with invariant SU(iV) elements. How the 
lattice Bianchi identities should be satisfied needs to be 
investigated [20l ]. 

The obvious difficulty with such an approach is the 
construction of operators — how should they appear in 
the plaquette representation? We conjecture that the 
building blocks of plaquette elements should be added 
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FIG. 9. Link correlation C(U, U') as a function of the number 
of updating sweeps at j3/d = 3 on a L D — 4 3 lattice for the 
MP □ and the HB Q algorithms. 
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FIG. 10. Plaquette correlation C(Q,Q') as a function of the 
number of updating sweeps at /3/d = 3 on a L D = 4 3 lattice 
for the MP □ and the HB O algorithms. 

together to form closed loops or bags of plaquettes, much 



as links are multiplied in their representation, so that all 
boundaries are closed (in the pure gauge theory). That 
construction is not unlike the depiction of bound cur- 
rents in the macroscopic treatment of electromagnetism, 
a topic which may be familiar to some readers [2l| . For 
the sectors of J PC with C = + for N > 2, the usual op- 
erator Re Tr Qc built from the trace of a product of links 
Qc = HcU(C) over a closed path C is obviously related 
to the Manton action for the element Qc, but for C = — 
a similar construction is not yet known to us. Investiga- 
tion of this approach is beyond the scope of the present 
work but remains a topic of interest. 



VII. CONCLUSIONS 

The availability of a general heatbath algorithm for 
SU(iV) lattice gauge theory has long been desired. While 
the Metropolis algorithm certainly is effective at thermal- 
ization, its implementation requires an overcommitment 
of resources to achieve a significant level of decorrelation. 
The advantage of the indirect heatbath algorithm is that 
a fewer number of trials are proposed for each link in a 
sweep, with consequent savings in time, energy, and ulti- 
mately money, as each flop exacts a price however small 
from the investigator. While we have not compared the 
indirect heatbath algorithm directly to that by Cabibbo 
and Marinari [3] , we expect to see a similar savings due to 
treating each element as a unit rather than decomposing 
it according to a covering set of subgroups. 

From the identification of the lattice action as the 
squared distance from the identity of the invariant pla- 
quette element we have constructed a general heatbath 
algorithm for SU(iV) pure gauge theory. The algorithm 
consists of proposing a number of new plaquette ele- 
ments, whose mean yields the updated link after dividing 
out the projection of the surrounding staple elements. 
The internal energy distribution obtained by the heat- 
bath algorithm agrees with that produced by Metropo- 
lis updating, and the link and plaquette decorrelation 
rates compare favorably, as does the execution speed. 
The indirect heatbath algorithm incorporates the Man- 
ton action, to which the Wilson action is a leading order 
approximation, accounting for the metric within the spe- 
cial unitary group rather than the general linear group. 
By breaking the relationship between the link to be up- 
dated and its replacement, the indirect heatbath algo- 
rithm achieves the maximal decorrelation of gauge in- 
variant objects with a minimal amount of effort. 
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